library(multiplot)
library(ggplot2)
library(dplyr)

setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/")
file=read.csv("20201120/at.least.one.AShM.in.DC.add.BF.beta0.add.CCHC.csv",head=T)
filea=read.csv("20201112做汇总表/all.FDR.sig.at.least.one.add.direction.same.diff.csv",head=T)
filea$id=paste(filea$Chr,filea$Start,sep = ":")
filea1=filea[filea$FDR.sig>1,]

file$id=paste(file$Chr,file$Start,sep=":")
file1=file[file$pattern.not.rm.dupl.num.DC>1,]
file2=file1[file1$BF_in_DC>1,]
file3=file1[file1$BF_in_DC>10,]


#思路20210124
test=file[file$pattern.nodupl.num.DC==1&file$pattern.nodupl.num.CC==1&file$pattern.nodupl.num.HC==1,]
test[test$pattern.nodupl.CC=="normal_down_tumor_down",]$pattern.nodupl.CC="normal_down-tumor_down"
test[test$pattern.nodupl.CC=="normal_up_tumor_up",]$pattern.nodupl.CC="normal_up-tumor_up"

test[test$pattern.nodupl.HC=="normal_down_tumor_down",]$pattern.nodupl.HC="normal_down-tumor_down"
test[test$pattern.nodupl.HC=="normal_up_tumor_up",]$pattern.nodupl.HC="normal_up-tumor_up"
test1=test[test$pattern.nodupl.DC==test$pattern.nodupl.CC&test$pattern.nodupl.DC==test$pattern.nodupl.HC,]

#原来的思路是从单发样本里8个pattern取交集，因此这22个位点筛选如下：
test2=file[file$DC.FDR.sig==12,]
test2=test2[test2$pattern.nodupl.num.DC==1,]

test3=file[file$CC.FDR.sig==8,]
test3=test3[test3$pattern.nodupl.num.CC==1,]

test4=file[file$HC.FDR.sig==8,]
test4=test4[test4$pattern.nodupl.num.HC==1,]

id=test2$unitID
file=file[file$unitID %in% id,]

normaldata=select(file,unitID,DC.con.beta0,DC.con.beta0.upper,DC.con.beta0.lower,DC.con.FDR,avsnp150)
normaldata$group=ifelse(normaldata$DC.con.beta0<0,"down","up")
normaldata.down=normaldata[normaldata$group=="down",]
normaldata.up=normaldata[normaldata$group=="up",]
normaldata.down=arrange(normaldata.down,DC.con.beta0)
normaldata.up=arrange(normaldata.up,DC.con.beta0)
normaldata=rbind(normaldata.up[normaldata.up$DC.con.beta0<0,],normaldata.down,normaldata.up[normaldata.up$DC.con.beta0>0,])

normaldata$num=1:dim(normaldata)[1]
normaldata$na=""
orders=select(normaldata,unitID,num)

p1=ggplot(normaldata,aes(y=DC.con.beta0,x=num))+scale_x_continuous(breaks = normaldata$num,labels = normaldata$na,position="top")+
  geom_point(aes(color=group),size=1.5)+geom_errorbar(aes(ymin=DC.con.beta0.lower,ymax=DC.con.beta0.upper,color=group),width=0.5)+labs(x="",y="con beta0")+
  geom_hline(yintercept = 0,color="red",linetype="dashed")+scale_color_manual(values = c("#3C5488","#E64B35"))+coord_flip()+theme_classic(base_size = 15)+
  theme(panel.grid.minor = element_blank(),panel.border = element_blank(),axis.line =element_line(size = 0.5))+theme(legend.position="none")

tumordata=select(file,unitID,DC.case.beta0,DC.case.beta0.upper,DC.case.beta0.lower,DC.case.FDR,avsnp150)
tumordata$group=ifelse(tumordata$DC.case.beta0<0,"down","up")
tumordata=merge(tumordata,orders,by="unitID")
tumordata=arrange(tumordata,num)
tumordata$na=""

p2=ggplot(tumordata,aes(y=DC.case.beta0,x=num))+scale_x_continuous(breaks = tumordata$num,labels = tumordata$na)+
  geom_point(aes(color=group),size=1.5)+geom_errorbar(aes(ymin=DC.case.beta0.lower,ymax=DC.case.beta0.upper,color=group),width=0.5)+
  labs(x="",y="case beta0")+geom_hline(yintercept = 0,color="red",linetype="dashed")+scale_color_manual(values = c("#3C5488","#E64B35"))+theme_classic(base_size = 15)+
  theme(panel.grid.minor = element_blank(),panel.border = element_blank(),axis.line = element_line(size = 0.5))+coord_flip()+theme(legend.position="none")


HC.case=read.table("./20201120/HC.case.beta0.txt",head=T,sep="\t")
HC.con=read.table("./20201120/HC.con.beta0.txt",head=T,sep="\t")

HC.case=merge(HC.case,orders,by="unitID")
HC.con=merge(HC.con,orders,by="unitID")
HC.case$na=""
HC.con$na=""

HC.case$group=ifelse(HC.case$HC.case.beta0<0,"down","up")
HC.con$group=ifelse(HC.con$HC.con.beta0<0,"down","up")

p4=ggplot(HC.con,aes(y=HC.con.beta0,x=num))+scale_x_continuous(breaks = HC.con$num,labels = HC.con$na,position="top")+
  geom_point(aes(color=group),size=1.5)+geom_errorbar(aes(ymin=HC.con.beta0.lower,ymax=HC.con.beta0.upper,color=group),width=0.5)+labs(x="",y="HC twin1 beta0")+
  geom_hline(yintercept = 0,color="red",linetype="dashed")+scale_color_manual(values = c("#3C5488","#E64B35"))+coord_flip()+theme_classic(base_size = 15)+
  theme(panel.grid.minor = element_blank(),panel.border = element_blank(),axis.line =element_line(size = 0.5))+theme(legend.position="none")

p5=ggplot(HC.case,aes(y=HC.case.beta0,x=num))+scale_x_continuous(breaks = HC.case$num,labels = HC.case$na)+
  geom_point(aes(color=group),size=1.5)+geom_errorbar(aes(ymin=HC.case.beta0.lower,ymax=HC.case.beta0.upper,color=group),width=0.5)+
  labs(x="",y="HC twin2 beta0")+geom_hline(yintercept = 0,color="red",linetype="dashed")+scale_color_manual(values = c("#3C5488","#E64B35"))+theme_classic(base_size = 15)+
  theme(panel.grid.minor = element_blank(),panel.border = element_blank(),axis.line = element_line(size = 0.5))+coord_flip()+theme(legend.position="none")

CC.case=read.table("./20201120/CC.case.beta0.txt",head=T,sep="\t")
CC.con=read.table("./20201120/CC.con.beta0.txt",head=T,sep="\t")

CC.case=merge(CC.case,orders,by="unitID")
CC.con=merge(CC.con,orders,by="unitID")
CC.case$na=""
CC.con$na=""

CC.case$group=ifelse(CC.case$CC.case.beta0<0,"down","up")
CC.con$group=ifelse(CC.con$CC.con.beta0<0,"down","up")

p7=ggplot(CC.con,aes(y=CC.con.beta0,x=num))+scale_x_continuous(breaks = CC.con$num,labels = CC.con$na,position="top")+
  geom_point(aes(color=group),size=1.5)+geom_errorbar(aes(ymin=CC.con.beta0.lower,ymax=CC.con.beta0.upper,color=group),width=0.5)+labs(x="",y="CC twin1 beta0")+
  geom_hline(yintercept = 0,color="red",linetype="dashed")+scale_color_manual(values = c('#708090','red'))+coord_flip()+theme_classic(base_size = 15)+
  theme(panel.grid.minor = element_blank(),panel.border = element_blank(),axis.line =element_line(size = 0.5))+theme(legend.position="none")

p8=ggplot(CC.case,aes(y=CC.case.beta0,x=num))+scale_x_continuous(breaks = CC.case$num,labels = CC.case$na)+
  geom_point(aes(color=group),size=1.5)+geom_errorbar(aes(ymin=CC.case.beta0.lower,ymax=CC.case.beta0.upper,color=group),width=0.5)+
  labs(x="",y="CC twin2 beta0")+geom_hline(yintercept = 0,color="red",linetype="dashed")+scale_color_manual(values = c('#708090','red'))+theme_classic(base_size = 15)+
  theme(panel.grid.minor = element_blank(),panel.border = element_blank(),axis.line = element_line(size = 0.5))+coord_flip()+theme(legend.position="none")

layout <- matrix(c(1,2,3,4,5,6), nrow = 1)
multiplot(plotlist=list(p1,p2,p7,p8,p4,p5),layout=layout)
